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1 .  INTRODUCTION 

The  physical  problem  we  are  addressing  in  this  study  is  the  rapid  deposi¬ 
tion  of  thermal  energy  on  the  surface  of  a  composite  plate  within  which  a 
delamination  exists.  The  high  energy  flux  will  stress  and  perhaps  vaporize  the 
surface  plies,  and  a  shock  stress  wave  is  initiated  through  thermal  expansion 
and  the  momentum  imparted  by  the  blowoff.  As  this  compressive  wave  initially 
passes  the  delamination,  the  crack  will  tend  to  close  and  the  wave  will  simply 
pass  through.  But  as  the  wave  reflects  off  the  bottom  free  surface,  it  will 
return  as  a  tensile  wave.  This  wave  will  open  the  delamination  and  produce 
excessive  stresses  at  the  crack  tip.  The  high  stress  field  could  conceivably 
cause  the  delamination  to  grow  and  result  in  the  catastrophic  failure  of  the 
plate.  In  addition  to  the  initiation  of  a  longitudinal  stress  wave,  a  shear 
wave  may  ensue  from  a  nonuniform  spatial  distribution  or  edge  effects  of  the 
laser  energy  flux.  Matters  are  further  complicated  by  the  anisotropic  and 
nonhomogeneous  makeup  of  the  composite  plate.  Hence,  the  longitudinal  and 
shear  waves  will  propagate  at  speeds  dependent  on  the  ply  orientation. 

Several  time  scales  may  be  identified.  The  energy  deposition  typically 
occurs  within  hundredths  of  a  microsecond.  Depending  on  the  wavelength  of  the 
source,  the  energy  may  be  deposited  on  just  the  surface  or  throughout  the 
plate,  but  in  either  case  it  results  in  an  almost  instantaneous  change  in  the 
temperature.  This  rapid  rise  in  temperature  may  produce  a  phase  change  at  the 
surface.  Subsequently,  the  pressure  in  the  gas  phase  becomes  very  high,  and 
the  adjacent  solid  portion  responds  with  a  stress  wave.  This  wave  will 
traverse  the  thickness  of  the  plate  on  the  order  microseconds.  The  conduction 
of  heat  takes  place  over  a  much  longer  time  scale.  For  example,  a  temperature 
change  of  1%  requires  on  the  order  of  100  milliseconds. 

1.1  The  Energy  Deposition  Phase 

The  analysis  tools  required  for  the  study  of  composite  structures  subjected 
to  vapid  thermal  pulses  are  conveniently  considered  in  three  parts,  two  of 
which  are  analyzed  in  some  detail  in  this  study.  The  first  part  is  the  deter¬ 
mination  of  the  nature  of  the  thermal  pulse  from  knowledge  or  assumptions 
about  the  source  of  that  pulse,  which  is  only  discussed  in  generalities  here. 
Sources  of  interest  include  laser  weapons,  nuclear  weapon  X-rays,  particle 
beams,  and  other  energy  sources.  The  analysis  of  this  part  of  the  problem 
must  consider  the  interaction  of  an  electromagnetic  source  with  the  material 
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j  of  the  composite  structure  and  must  analyze  the  subsequent  radiation  transport 

in  that  material.  Analysis  techniques  are  available,  but  experimental  data 
sufficient  to  define  material  properties  for  energy  deposition  in  polymeric 
materials  characteristic  of  composites  are  sparse  or  nonexistent. 

The  outcome  of  a  study  of  this  radiation  transport  is  a  time-position 
profile  of  the  thermodynamic  state  of  the  composite  structure.  Fortunately, 
to  study  and  develop  the  tools  required  for  the  remaining  two  parts  of  the 
problem,  it  is  not  necessary  to  have  specific  numerical  results.  Rather,  it 
is  only  necessary  to  know  the  general  nature  of  the  results  so  that  the  tools 
!  for  the  rest  of  the  analysis  are  sufficient  to  cope  with  the  range  of  possibi¬ 

lities.  A  short  discussion  of  the  general  nature,  and  of  the  differences  due 
to  the  variety  of  possible  source  types,  is  given. 

It  is  convenient  to  classify  energy  sources  according  to  their  range  of 
initial  spatial  influence  and  their  temporal  duration.  The  best  studied 
source  is  the  detonation  of  a  nuclear  weapon  near  a  structure  of  interest. 

The  predominant  energy  of  a  modern  fusion  device  is  in  a  substantial  X-ray 
output  that  impinges  on,  and  is  absorbed  into,  the  adjacent  material.  If  the 
detonation  is  in  the  atmosphere,  the  surrounding  air  absorbs  the  energy,  and  a 
blast  wave  is  formed  that  propagates  outward  and  strikes  nearby  structures. 
Alternatively  (and  exclusively  if  there  is  no  significant  atmosphere  present), 
the  radiation  directly  reaches  a  nearby  structure.  The  energy  is  then  absorbed 
by  the  material  in  the  structure. 

The  time  scales  of  this  energy  transport  phase  are  typically  very  short  in 
comparison  to  the  other  significant  time  scales  of  the  problem.  In  many  cases, 
it  is  considered  instantaneous.  The  characteristic  depth  of  the  energy  deposi¬ 
tion  depends  markedly  on  the  composition  of  the  structure.  Materials  with  high 
atomic  numbers  absorb  the  energy  in  very  shallow  depths,  while  so-called 
"low-Z"  materials  have  much  greater  absorption  depths. 

Thus,  energy  deposition  due  to  nuclear  devices  is  typically  of  very  short 
duration.  The  depths  of  deposition  can  be  either  thin  compared  to  structural 
dimensions  (e.g,,  a  skin  thickness  of  a  missile)  or  of  that  same  size. 

When  considering  laser  weapons,  there  are  important  differences  in  this 
picture  due  to  the  longer  wavelengths  of  the  radiation.  The  time  scales,  for 
a  single  pulse,  are  still  short.  However,  the  deposition  thickness  is  very 
small  and,  as  a  result  of  the  very  high  energy  densities,  it  will  vaporize  and 
"blow  off"  on  a  short  time  scale,  which  can  have  a  significant  quenching  effect 
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on  the  energy  deposition  itself.  In  this  case,  there  can  be  important  inter¬ 
actions  between  the  thermodynamics  of  the  material,  the  resulting  wave  actions, 
and  the  energy  deposition  itself.  For  that  reason,  multiple  and  rapidly  pulsed 
weapon  outputs  are  under  consideration.  The  tools  to  be  developed  here  must 
be  able  to  handle  this  important  case. 

1 . 2  The  Initial  Thermodynamic  State 

The  result  of  an  analysis  of  the  energy  deposition  and  radiation  transport 
is  a  specified  physical  state  for  all  material  points  in  the  structure  as  a 
function  of  time.  As  was  discussed,  the  time  scale  of  this  phase  of  the 
problem  is  typically  very  short  and,  for  the  present,  will  be  considered  to  be 
instantaneous . 

The  "physical  state"  means  the  thermodynamic  state.  The  nature  of  this 
state  can  be  described  by  referral  to  a  typical  phase-state  diagram  for  a 
material.  Figure  1  shows  a  typical  temperature-density  plot  of  all  equilibrium 
thermodynamic  states  for  a  metal.  The  boundaries  between  the  solid,  liquid, 
and  vapor  phases  are  shown.  Also  shown  is  the  locus,  of  points  at  a  constant 
one-atmosphere  pressure.  Along  that  locus,  the  melt  and  vapor  points  are  iden¬ 
tified.  The  critical  point  and  triple  line  are  also  indicated. 

Only  certain  parts  of  this  diagram  are  of  importance  to  the  present 
analysis.  A  typical  material  point  is  initially  at  standard  temperature  and 
pressure.  For  clarity,  Figure  2  will  be  used  for  this  discussion;  this 
initial  point  is  labeled  as  point  A. 

For  instantaneous  energy  deposition,  the  temperature  (and  also  the  internal 
energy)  is  suddenly  increased.  Since  there  is  insufficient  time  for  material 
motions,  the  material  will  remain  at  constant  mass  density.  Consequently,  the 
material  point  will  achieve  the  state  at  point  B,  a  point  directly  above  point 
A.  Depending  on  the  magnitude  of  the  energy  absorbed,  this  point  can  be  well 
into  the  liquid  or  vapor  regions.  However,  since  the  mass  density  remains  at 
the  initial  value,  the  pressure  will  be  very  high.  As  an  example,  point  B  is 
shown  along  a  one-megabar  pressure  line,  corresponding  to  a  point  typically 
well  above  the  critical  point. 

We  will  discuss  in  detail  the  subsequent  material  motions  and  wave  propaga¬ 
tions  that  result  from  such  initial  states  throughout  a  structure.  A  material 
point  will  not  remain  at  poinl  B,  but  the  material  will  expand  and  the  pressure 


will  reduce  eventually  back  to  normal  pressure.  Since  these  motions  occur 
after  the  energy  deposition,  they  are  adiabatic  and,  as  a  consequence,  the 
path  followed  in  this  thermodynamic  phase  space  is  along  an  isentrope.  A 
typical  isentrope  from  point  B  is  indicated.  The  one  shown  happens  to 
intersect  the  vapor-liquid  dome,  but  it  will  not  if  point  B  is  at  a 
sufficiently  high  temperature.  The  final  states  of  the  material  will 
ultimately  be  at  very  low  densities,  and  the  material  will  be  hot  and  expanded. 

This  discussion  serves  to  identify  the  thermodynamic  description  and  model 
that  is  needed  for  the  initial  deposition  phase  for  an  analysis  of  the  type 
studied  here.  A  relatively  precise  description  of  all  states  at  nominal  and 
lower  densities  is  needed.  A  model  of  the  phase  boundaries  and  the  shape  of 
the  pressure  curves  and  adiabats  (isentropes)  in  this  region  is  also  needed. 

Subsequent  wave  motions  can  introduce  further  states.  In  the  case  that  a 
layer  of  surface  material  vaporizes  and  blows  off  at  high  velocity,  or  a  thin 
layer  of  material  is  spalled  off,  then  a  shock  wave  will  be  generated  that 
propagates  into  the  interior  of  the  structure.  This  shock  wave  is  similar  to 
one  that  would  be  generated  by  an  impact  at  the  surface.  The  resulting  shocked 
states  lie  along  a  Hugoniot  curve  centered  at  the  initial  point  of  the 
material.  A  typical  Hugoniot  curve  from  the  standard  temperature  and  pressure 
is  also  shown  in  Figure  2.  We  can  see  that  these  states  are  at  higher  than 
normal  density  and  temperature. 

In  the  case  of  composite  materials,  specific  data  on  these  thermodynamic 
states  are  very  sparse.  In  the  present  Phase  I  study,  a  thorough  literature 
search  has  not  been  conducted,  but  preliminary  investigations  uncovered  only  a 
very  few  thermodynamic  properties,  including  only  the  initial  density,  the 
Gruneisen  parameter,  the  wave  speed,  and  thermal  expansion  at  standard 
temperature  and  pressure.  No  data  were  found  on  temperature  states  above  a 
few  hundred  degrees,  and  none  describing  phase-change  mechanisms.  Clearly, 
before  we  can  have  sufficient  confidence  in  energy  deposition  studies,  much 
more  work  is  needed  in  this  area:  both  analytical  descriptions  of  the 
thermodynamics  of  composite  materials  and  experimental  tests  to  verify  and 
calibrate  those  analytical  models.  The  additional  features  arising  from  the 
multiple-constituency  of  the  composite  materials  will  also  require  further 
investigation. 
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sical  Wave  Profiles 


To  understand  the  computer  results  to  be  shown  shortly,  we  first  describe 
the  general  nature  of  the  stress  pulses  that  will  arise  in  the  cases  of  rapid 
thermal  pulses.  An  idealized  description  is  given  here,  followed  by  a 
presentation  of  some  actual  code  output  that  includes  all  of  the  real  physics 


of  the  problem. 

For  an  idealized  case,  consider  a 
plate  of  thickness  t  in  the  x  direction. 
We  may  assume  that  there  is  an 
instantaneous  triangular  uniform  energy 
deposition  at  the  free  surface  of  that 
plate,  with  the  maximum  value  at  the 
surface  and  dropping  linearly  to  zero 
at  some  depth  h  less  than  t. 

There  is  a  resulting  high  pressure 
in  that  deposition  region  that  produces 
an  initial  pressure  profile  as  shown  in 
sketch  1 . 

The  general  effect  of  this  initial 
state  can  be  understood  by  considering 
the  special  case  of  a  linear  elastic 


material  for  which  superposition  holds . 


In  that  case,  the  bal;  ice  of  linear 
momentum  reduces  to  the  well-known  one¬ 
dimensional  wave  equation.  If  the 
pulse  is  not  near  the  free  surface, 
then  the  initial  pulse  shown  above 
would  split  into  two  equal  parts,  one 
traveling  to  the  right  at  the  wave 
velocity  c,  and  one  to  the  left.  Thus, 
a  short  instant  At  later,  these  two 
waves  look  as  shown  in  sketch  2. 

However,  there  is  a  free  surface  at 
x  =  0  to  consider.  The  condition  that 
the  pressure  be  zero  at  x  -  0  is  satis¬ 
fied  by  adding  a  third  wave,  opposite  to 


the  part  traveling  left  in  shape  and 
sign  (sketch  3),  that  will  at  all 
times  cancel  that  compressive  puls®. 
The  result  of  these  three  waves  at 
the  time  At  is  shown  in  sketch  4. 

After  the  pulses  all  clear  the 
free  surface,  the  result  is  a 
classical  "N-shaped"  wave  as  shown 
in  sketch  5.  Therefore,  in  this 
simplistic  analysis,  we  can  see 
that  the  net  effect  of  an  initial 
triangular  energy  deposition  is  a 
wave  with  a  triangular  compressive 
front  followed  by  a  triangular 
tensile  tail  of  equal  magnitude. 
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This  analysis  assumes  instantaneous  energy  deposition.  A  slightly  modified 
picture  occurs  when  the  deposition  time  t  is  not  small  compared  to  the  pulse 
depth  h  divided  by  the  wave  speed  c.  In  that  case,  for  the  same  total  energy 
and  deposition  depth,  the  final  wave  after  deposition  time  will  have  a  reduced 
peak  and  will  have  a  total  spatial  width  of  h  +  ctQ.  However,  there  will  still 
be  a  following  tensile  tail  of  the  same  magnitude  as  the  compressive  front. 

Real  materials  are  not  linear  elastic,  particularly  at  the  high  temperature 
and  pressure  states  encountered  during  the  energy  deposition  phase.  The  most 
important  feature  of  nonelastic  actual  behavior  wili  be  a  finite  tensile 
strength,  due  either  to  the  inherent  strength  of  the  solid,  or  due  to  the 
reduced  strength  of  a  vaporized  or  partially  vaporized  state. 

Assume,  for  example,  that  the  energy  deposition  magnitude  is  such  that  the 
peak  compressive  stress  in  the  above  sketches  is  5  kilobars,  and  that  the 
tensile  spall  strength  is  1  kilobar.  Then,  at  some  time  between  the  first 
growth  of  the  tensile  tail  at  the  free  surface  and  the  time  when  the  N-shaped 
wave  would  have  left  the  free  surface,  the  tensile  stress  at  some  distance  in 
from  the  free  surface  will  reach  1  kilobar,  and  tensile  spall  will  occur. 
Indeed,  with  the  values  stated  in  this  simple  example,  that  initial  spall 
thickness  will  be  exactly  one-fifth  of  the  deposition  depth  h.  That  reduces 
the  stress  at  the  spall  plane  back  to  zero,  and  a  tensile  pulse  will  again 
begin  to  grow,  resulting  in  a  second  spall  of  equal  thickness.  In  the  final 
analysis,  the  resulting  wave  propagating  into  the  structure  will  have  a  tensile 
tail,  but  limited  in  magnitude  to  the  tensile  spall  strength.  In  the  example 
here,  that  tail  would  be  limited  to  1  kilobar. 

The  final  stress  wave  profile  will,  of  course,  depend  on  the  exact  nature 
of  the  energy  deposition.  In  a  one-dimensional  thermoelastic  analysis  neglect¬ 
ing  vaporization  and  tensile  strength,  Paramasivam  and  Reismann  (1986) 
(Reference  7)  find  a  similar  compressive/ tensile  wave  resulting  from  a 
Gaussian  deposition  in  space  and  time. 

1.4  Assumptions 

A  complete  analysis  of  this  complicated  thermomechanical  problem  is  beyond 
the  scope  of  this  Phase  I  work,  but  there  are  a  number  of  assumptions  that  may 
be  made  for  the  problem  to  become  tractable  for  analysis.  First,  we  tacitly 
assume  that  the  energy  flux  is  large  enough  to  vaporize  the  first  few  plies 
but  is  not  sufficiently  strong  to  vaporize  to  a  significant  depth  in  the 
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plate.  If  we  assume  that  the  delamination  does  not  reside  too  close  to  the 
upper  surface,  then  the  effects  of  the  delamination  will  not  be  strongly 
coupled  to  the  shock  wave  initiation.  Hence,  the  interaction  of  the  stress 
wave  with  the  delamination  may  be  treated  separately  from  the  physics  of  its 
initial  generation.  Furthermore,  if  we  wish  to  analyze  the  stresses  at  the 
delamination  in  a  planar  setting,  we  are  forced  to  make  some  assumptions 
regarding  the  areal  extent  of  the  energy  deposition  and  the  geometry  of  the 
delamination.  We  may  assume  either  that  the  energy  is  uniformly  distributed 
over  an  area  much  larger  than  the  characteristic  width  of  the  delamination,  or 
that  the  energy  is  deposited  uniformly  over  a  thin  but  very  elongated  region. 
The  first  assumption  seems  more  realistic,  but  the  second  allows  us  to  study 
effects  related  to  the  position  of  the  energy  deposition  with  respect  to  the 
delamination.  The  delamination  itself  must  be  regarded  as  having  some  large 
extension  down  the  plate.  In  other  words,  the  delamination  is  a  long  cavity 
with  small  width  and  even  smaller  thickness.  With  these  assumptions,  the 
deformation  is  uniform  along  the  long  axis  of  the  delamination  (and  the  energy 
deposition  area),  and  the  plate  exists  in  a  state  of  plane  strain  (Figure  3). 

The  assumption  of  a  large  uniform  energy  deposition  further  reduces  the 
dynamics  of  the  stress  wave  propagation  to  a  one-dimensional  problem,  since  the 
effects  of  the  crack  do  not  interact.  Also,  since  the  uniaxial  stress-strain 
relation  is  independent  of  the  ply  orientation,  the  stress  wave  propagates  as 
in  a  homogeneous  material.  It  is  only  after  the  stress  wave  interacts  with 
the  delamination  that  the  orthotropic,  nonhomogeneous  properties  become 
apparent.  Accordingly,  if  we  assume  that  the  initial  stress  wave  is  entirely 
compressive,  and  the  closed  delamination  cannot  slip  (infinite  friction),  then 
the  compressive  wave  will  pass  through  unchanged. 

With  these  assumptions,  the  problem  may  be  approached  with  the  two  existing 
codes,  FEAPICC,  a  two-dimensional,  finite  element  code  for  propagating  cracks 
between  differing  anisotropic  materials,  and  WONDY,  a  one-dimensional,  finite 
difference  code  modeling  the  thermomechanics  of  the  energy  deposition.  Before 
moving  on  to  a  discussion  of  the  simulations  using  these  codes,  the  next 
section  presents  more  details  of  the  theoretical  model. 
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2.  TKEC^STICAL  DEVELOPMENTS 
2.1  Continuum  Mechanics  Overview 

We  present  here  a  brief  description  of  the  continuum  model  with  emphasis 
on  particular  aspects  that  relate  to  the  deposition  of  energy  in  composite 
laminates . 

The  motions  are  governed  by  the  general  equations  of  continuum  mechanics, 
as  are  used  in  all  studies  of  fluid,  solid  and  structural  mechanics.  The 
computer  codes  used  in  this  and  other  studies  transcribe  those  equations  into 
an  approximate  form  suitable  for  solution  on  a  computer. 

The  equations  are  recorded  here  in  Lagrangian  form.  The  position  of  a 
material  point  in  a  reference  Lagrangian  coordinate  system  is  denoted  by  x°. 

The  spatial  position  x  is  then  a  function  of  x°  and  time  t: 

X  *  X  (x°,  t)  .  (1) 


A  variety  of  kinematic  definitions  include  those  for 


Velocity: 

Acceleration: 

Deformation  gradient  tensor: 
Velocity  gradient  tensor: 
Stretching  tensor: 


v  »  x 

a  «  v 


F  *  V°x 

***  *»*  M 

#  “l 

L  -  F  F 


2D  «  L  +  L 


(2) 


Here  the  superposed  dot  denotes  the  material  time  derivative  and  V°  is  the  del 
operator  with  respect  to  x°.  The  stress  tensor  is  denoted  by  g,  which,  to  be 
precise,  is  the  first  Piola-Ki’chhof f  tensor  related  to  the  more  standard 
Cauchy  stress  tensor  T  by  the  relation 


with 


-1  T 

g  -  J  T  (f  V 

J  =  det  F 


The  pressure  p  is  given  by 


(3) 

(4) 


p  =  -1/3  tr  g 


Using  the  symbol  p  to  denote  the  mass  density  and  pQ  the  density  in  the 


(5) 


reference  configuration,  the  balance  of  mass  can  be  written  as 

p  +  pV*v  =  0  ,  (6) 

the  balance  of  linear  momentum  as 

V°*g  +  pQ*b  =  pQ«a  ,  (7) 

and,  finally,  the  balance  of  energy  as 

pe  =  tr(TD)  -  V»<j  +  pr  ,  (8) 

where  b  is  the  body  force,  ^  is  the  heat  flux  vector  and  r  is  the  deposited 
energy  per  unit  mass. 

These  equations  are  supplemented  by  the  constitutive  equations  describing 
the  material.  There  is  a  relation  between  the  internal  energy  e,  the  density  p 
and  the  pressure  p,  which  we  write  as 

p  =  p(e ,  p)  (9) 

While  this  single  equation  of  state  form  is  sufficient  to  describe  completely 
the  thermodynamics,  it  is  more  common  to  also  include  descriptions  for  the 
temperature  T  and  entropy  q: 

T  -  T(e,  p) 

(10) 

0  ■  n(e .  p) 

Then  these  three  thermodynamic  equations  can  be  solved  to  use  any  two  thermo¬ 
dynamic  variables  as  independent.  The  discussions  of  the  previous  section  used 
T  and  p;  the  corresponding  forms  are  given  by 

p  =  p(p,  T) 

e  =  e  (  p,  T)  (11) 

n  =  n(p»  t) 

and  constant  pressure  or  constant  entropy  states  lie  along  some  curve  in  p  -  T 
space  as  described  above.  Phase  boundaries  become  curves  in  this  p  -  T  space. 

To  complete  the  description,  a  relation  for  the  shear  stress  components  is 
needed.  This  usually  takes  the  form  of  equations  for  the  stress  deviator  tensor: 


oD  =  o  -  p  fi 


Relating  the  rate  of  change  of  to  the  stretching  tensor  D  via  a  shear 
modulus  G  is  usual.  The  shear  modulus  G  can  depend  on  the  thermodynamic  state, 
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i.e.,  the  temperature  and  density.  Plasticity  relations  are  used  to  limit  the 
values  of  the  stress  deviators  when  plastic  flow  occurs;  the  yield  strength  Y  is 
also  dependent  on  the  thermodynamic  state.  Polymeric  materials  of  composite 
structures  may  require  more  general  and  time-dependent  models.  Again,  little  is 
presently  known  about  suitable  models.  The  most  that  appear  to  be  available  are 
quasi-t  ;atic  dependences  of  strength  and  stiffness  versus  temperature,  and  only 
to  the  few  hundred  degrees  for  which  composite  materials  maintain  their  primrry 
structural  integrity.  Finally,  fracture  criteria  are  needed  to  describe  spall 
(tensile  failures)  and,  in  the  case  of  composite  structures,  ply  delaminations. 

It  is  clear  that  there  is  an  imposing  amount  of  information  that  goes  into 
the  final  solution  of  these  problems.  Fortunately,  there  are  well  developed 
and  tested  finite  difference  codes,  such  as  the  one-dimensional  code  WONDY  and 
the  two-dimensional  code  CSQ  that  have  been  used  in  the  present  program. 
Provisions  for  a  number  of  different  material  model  types  are  available,  as 
well  as  both  analytical  and  tabular  descriptions  of  the  thermodynamics, 
including  complete  three-phase  boundary  descriptions  and  transitions.  At  the 
present,  for  normal  wave-propagation  studies,  it  must  be  said  that  the  models 
are  better  developed  than  warranted;  however,  the  experimental  data  to  provide 
inputs  to  these  models  (or  even  to  simply  discover  what  type  of  a  model  is 
appropriate)  are  lacking.  Furthermore,  even  for  generic  input  models,  the 
effects  of  composite  laminations  or  delaminations  on  the  wave  propagation 
mechanics  has  not  been  considered.  It  is,  of  course,  that  issue  that  is  the 
primary  focus  of  the  present  study. 

2.2  Singular  Finite  Element  Formulation 

In  the  present  study,  singular  finite  elements  are  used  near  the  tip  of  a 
delamination  crack  to  account  for  stress  singularity  at  the  crack  tip.  These 
elements  incorporate  stress  and  displacement  fields  from  the  closed-form 
solutions  and  therefore  are  extremely  accurate  and  efficient.  The  singular 
elements  are  then  combined  with  the  regular  isoparametric  elements  in  the 
surrounding  region  so  that  the  standard  finite  element  procedures  can  be  used 
to  obtain  displacement  at  each  node. 

The  development  of  singular  elements  is  based  on  a  hybrid  functional 
derived  in  Tong  et  al.  (1973)  (Reference  10).  This  hybrid  functional  was  used 
by  Lin  and  Mar  (1976)  (Reference  5)  for  the  study  of  bi-material  crack  problems 
and  recently  by  Aminpour  and  Holsapple  (Aminpour,  1986)  (Reference  1)  for  the 
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dynamic  analysis  of  cracks  between  two  anisotropic  materials.  For  a  two- 
dimensional  continuum  divided  into  m  elements,  the  hybrid  functional  can  be 


defined  as 


where 


n  *  l  n 


=  J  |  f  e*(g  +  g°)  -  pu  *b  + 

Y  v"  f  -  \ 

+  /  (u  “  VJ  *X  ds  ~  /  u*t  ds  > 

Js.  ~  JS  ~  ~  \ 

m  o  » 


,  1  •  •  I 

-  pu  *b  +  -  pu  *u  J 


In  the  above,  a,  e,  a  ,  and  u  are  the  stress,  strain,  thermal  stress  and  dis¬ 
placement,  respectively.  Otner  variables  are  the  body  force  b,  the  density  p 
and  the  tractions  x  and  t  on  the  boundaries.  The  displacement  v  along  the 
interelement  boundary  is  assumed  independently  of  the  interior  displacement 
u.  The  constraint  integral  over  is  added  to  enforce  the  continuity  of  the 
displacements  between  the  singular  element  and  a  regular  element  that  uses 
different  interpolation  functions  for  the  displacement  components. 

Making  use  of  the  stress,  strain,  and  displacement  fields  obtained  from  the 
closed-f  m  solution  of  a  semi-infinite  crack,  we  can  assume 

u  =  UB 

a  =  AEIJ  +0°  (15) 


T  =  Rft 

where  the  (S's  are  unknown  constants  to  be  determined  from  the  finite  element 
solution  of  the  overall  problem. 

If  the  field  variables  in  Equation  (15)  are  used,  then  all  elasticity 

equations  are  satisfied  exactly  at  any  time  t,  and  the  volume  integral  in  the 

hybrid  functional  can  be  reduced  to  a  boundary  integral.  The  Euler  equation 

for  the  variational  functional  is  simply 

u  =  v  on  S. 

~  m 

and  (16) 


X  =  n  »o 

where  n  is  the  normal  vector. 


on  S 

o 


The  interelement  boundary  d is  pi a cement  u  can  be  assumed  in  the  following  form: 


wrja-jij'^xwTjcrjn^wTuciJtj.  .vwjww  xn  *.-  x.-.  k 


in  which  c^  is  the  nodal  displacement  vector,  and  L  is  chosen  such  that  u  is 
commit,  '-le  with  the  surrounding  regular  elements.  For  example,  L  may  be 
chosen  to  vary  linearly  between  two  nodes  along  the  interelement  boundary  if 
constant  strain  elements  are  used  in  the  surrounding. 

Substituting  the  assumed  quantities  from  Equations  (15)  and  (17)  into 
the  hybrid  functional  and  taking  the  variation  of  the  functional  with  respect 
to  (3,  we  obtain  (Aminpour,  1986)  (Reference  1) 


~1. 


(3  =  B  q  with  B  =  P  G 


and 


where 


/lT  1  »T  •  »T  T  \ 

^  a  H  -  2  a  "a  -  a  va  -  a  f  j 


K  =  BTHB  -  BTM.B  -  BTM,B  -  2BTM-B 
~  ~  ~~  ~  ~1~  ~  ~2~  ~  ~3~ 


dt 


M  =  B  M2B 
V  =  btm  b  +  btm.b 

~  ~  ~2  ~  ~  ~  i  ~ 

T 

F  =  B  F 
~  ~  ~s 


Ti: a  remaining  matrices  are  defined  as  follows: 


ETAE  dv 


H-2  "k 

P-'/s-  -T"- 


m 


m 


m 


m 


m 


(UTpb  -  Eg°)  dv 


•  T  • 

U  pu  dv 


UTpU  dv 


•T 

U  pU  dv 


T 

R  L  ds 


ds 
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(18) 

(19) 


(20) 


(21) 
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in  which  the  dots  denote  the  tine  derivatives  and  A  is  the  stress-strain 
relation  matrix.  The  matrix  G  is  defined  as 


E  * 


3_ 

3x 

0 

9__ 

9y 


0 

3_ 

3y 

9__ 

3x 


(22) 


In  Equation  (19),  K  is  the  element  stiffness  matrix,  M  is  the  mass  matrix,  V 

is  the  damping  matrix,  and  F  is  the  element  force  vector.  These  quantities 

can  be  assembled  with  regular  elements  using  standard  finite  element  procedures. 

The  summation  of  it  in  Equation  (19)  over  the  entire  domain  and  the  use  of 
m 

variational  theorems  will  yield  the  following  governing  differential  equation 
for  elasto-dynamic  problems: 

Mq  +  Vq  +  Kq  =  F  (23) 


where  M,  V,  K,  and  F  are  assembled  quantities.  Although  there  is  no  damping 
considered  in  this  problem,  the  matrix  V,  called  the  "pseudo-damping"  matrix,  is 
present.  Also,  for  a  propagating  crack,  the  matrix  M  is  symmetric,  while  the 
matrices  K  and  V  are  not,  as  can  be  seen  from  Equations  (20)  and  (21).  For  a 
stationary  crack,  the  matrix  V  vanishes  and  the  matrix  K  becomes  symmetric.  The 
I  existence  of  the  pseudo-damping  matrix  V  and  the  nonsymmetry  of  the  K  matrix 

;  occur  only  in  the  formulation  of  singular  elements.  These  complexities  arij.e 

[because  the  eigenfunctions  for  the  singular  element  were  derived  with  respect 
to  a  moving  local  coordinate  system  at  the  crack  tip;  therefore,  the  matrices 
M,  V,  and  K  are  functions  of  time. 

|  The  derivations  of  element  stiffness  and  mass  matrices  for  the  surrounding 

i 

regular  elements  can  be  found  elsewhere,  for  example,  in  Aminpour  (1986) 
(Reference  1). 

|  2.3  Strain  Energy  Release  Rate 

The  strain  energy  release  rate  is  defined  as  the  energy  released  per  unit 
of  new  crack  surface  generated  by  a  propagating  crack.  For  a  general  two-crack 
problem,  the  strain  energy  release  rate  can  be  expressed  by  the  following  crack 
closure  integral  (Lawrence  and  Masur,  1971)  (Reference  3): 


-  H»i  f 

6*0  5  A 


0  (x,0)[w(x-6,0)] 
z 


T  (x,0)[u(x-6,0)] 
xz 


dx 


(24) 
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where  the  first  term  in  the  integral  can  be  defined  as  (mode  1)  and  the 
second  term  as  (mode  2).  The  stresses  are  evaluated  at  a  distance  x  ahead 
of  the  crack  tip  and  the  corresponding  displacements  are  calculated  at  a  distance 
6  -  x  behind  the  tip  of  a  crack. 

In  the  present  study,  the  stress  and  displacement  are  obtained  from  an  eigen¬ 
function  expansion  method  with  each  term  accurate  to  a  constant  B^.  The  B 
matrix  is  then  determined  by  matching  the  local  crack-tip  solution  (singular 
elements)  with  far-field  solutions  (regular  elements).  Once  B's  are  found,  the 
stresses  and  displacements  are  known,  and  Equation  (24)  can  be  integrated 
numerically.  Note  that  since  the  stress  or  strain  is  oscillating  near  the  tip 
of  an  interface  crack  between  two  different  materials,  the  and  results  as 
obtained  from  Equation  (24)  generally  depend  upon  the  6  values  chosen  (Froula 
et  al.,  1980)  (Reference  2).  Therefore,  the  usefulness  of  separate  G^,  G^ 
values  for  bi-material  crack  problems  remains  to  be  studied.  However,  the  total 
strain  energy  release  rate  G^  has  been  shown  to  be  independent  of  6  (Walkar 
and  Lin,  1987)  (Reference  11). 


2.4  Plane  Strain 

Consider  a  lamina  whose  fibers  lie  ir.  the  XY  plane  as  shown  in  Figure  3. 
Consider  this  ply  to  be  a  monoclinic  material  with  a  stress/strain  relation  of 
the  form 


(25) 


where  is  the  compliance  matrix  and  is  the  coefficient  of  thermal  expan¬ 
sion.  Here,  we  use  the  contracted  notation  for  stress  and  strain  and  T  is  the 


temperature  above  the  reference. 

Assume  that  the  deformations  are  independent  of  the  X  direction  so  that 


y  =0 
'xy 


X  'xz 

for  which  we  may  solve  for  stresses  acting  on  the  X  face. 


(26) 
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°x  S  d  (S16S26  *  S66S12)ay  +  (S16S36  "  S66S13)az  +  (S16°fc  '  S66al)T 
Txy  =  I  (S12S16  “  SllS26)oy  +  (S16S13  ’  S36Sll)az  +  ^S16°l  "  Slla6)T  (27) 


Txz  "  S55  Tyz 


d  =  S11S66  -  S16  *  (28) 

Substituting  these  expressions  back  into  the  three-dimensional  stress/strain 
relations  and  adopting  the  two-dimensional  stress/strain  definition  as 

(  S  \  fbn  bi2  0  1  (  °»  W  “y  \ 


b12  b22  0 


°z  \  *  1  az  V  • 


0  0  b, 


bll  "  S22  +  (2  S 12S 16S 26  “  S66S12  "  SllS26)/d 
b22  =  S33  +  (2  S13S16S36  "  S66S13  ”  SllS36)/d 


b  12  ~  S23  +  (S16S12S36  +  S16S26S13  "  S66S13S12  "  SllS26S36)/d 


b33  =  S44  "  S45/S55 


^S26S16  "  S66S12^  ( 

°v =  a2 +  \ - 3 - /  “1  +  V 

fS36S16  "  S66S13  ^  ( 

az  =  a3  *  V - d - )  «1  +  \ 


S12S16  “  S26S11 


S13S 16  _  S  3  6S 1 1 


Note  that  under  the  plane  strain  assumption,  the  material  behavior  becomes 
orthotropic.  This  is  as  expected,  since,  depending  on  the  orientation  of  the 
fibers,  the  YZ  plane  cuts  the  fibers  in  circles,  ellipses,  or  straight  lines, 
as  seen  in  Figure  4. 


xv  j*v  x\j  x  y  ^rvt  jgs  km  ^ru^aj^iLK^K  x  < 


[iH(tJUifliOSlLH*:flX«XnX?l*)\>lAX/lAJ..  aUUUUJ  MJ'j'.JWlAn JV1  ATJmJWJrA.* 


Figure  4.  Plane  Sections  of  a  Laminate 


The  inverse  relationship  is  given  as 


where 


*  b„  /a 

a.  ^ 

11 

22 

12 

22 

■  bn/a 

a33 

/  a 


b  b 
11  22 


(31) 


(32) 
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3.  FINITE  ELEMENT  VALIDATION  TEST 

In  order  to  carry  out  this  study,  the  FEAPICC  code  was  transported  from  a 
CDC  60-bit  computer  to  the  Apollo  32-bit  machine.  This  necessitated  converting 
the  code  from  single  precision  to  double  precision.  However,  since  no  double¬ 
precision  complex  arithmetic  is  available  on  the  Apollo,  most  of  the  singular 
element  computations  were  left  in  complex  form  with  no  apparant  loss  in 
precision.  Conversion  to  double  precision  was  especially  important  for  dynamic 
problems  and  problems  with  many  unknowns. 

In  addition  to  these  conversion  problems,  new  theoretical  developments  and 
applications  and  new  I/O  formats  necessitated  other  modifications  that  required 
testing.  This  section  documents  tests  performed  to  validate  the  accuracy  of 
the  code. 


3.1  Singular  Element  Te^t 

This  section  investigates  the  accuracy  of  the  singular  element  when  its 
dimensions  become  small  compared  to  the  size  of  the  delamination.  A  typical 
ply  thickness  of  a  composite  is  on  the  order  of  0.005  inch,  while  the  delamina¬ 
tion  itself  may  be  as  large  as  1  inch.  Although  nondestructive  tests  would 
rule  out  such  a  large  delamination,  it  is  conceivable  that  a  delamination  could 
grow  to  such  lengths  under  normal  operating  conditions.  To  represent  the 
stress  variation  properly,  one  or  two  elements  would  be  required  to  lie  within 
each  ply.  We  find  that  the  ratio  of  the  crack  length  to  the  singular  element 
size  may  be  very  large.  A  simple  static  analysis  shows  that  the  stresses  grow 
by  the  square  root  of  the  ratio  of  the  crack  length  to  the  dis¬ 
tance  to  the  crack  tip.  If  the  singular  element  occupies  such  a  small  region 
around  the  crack  tip,  this  means  that  regular  isoparametric  elements  next  to 
the  singular  element  must  resolve  large  stress  gradients. 

O'Leary  (1981)  (Reference  6),  in  a  one-dimensional  error  analysis  and 
supporting  numerical  studies,  shows  that  as  the  singular  element  becomes  small, 
the  error,  in  the  energy  norm,  converges  to  zero.  However,  the  rate  of  conver¬ 
gence  suddenly  slows  after  the  singular  is  reduced  below  a  certain  size.  The 
rate  of  convergence  then  becomes  that  which  would  be  expected  if  no  singular 
element  was  employed.  The  error  in  the  stress  intensity  factor  shows  similar 
behavior  except  that  no  convergence  is  observed  after  the  singular  element  is 
reduced  below  a  certain  size.  In  any  case,  the  errors  were  significantly 
improved  by  the  use  of  the  singular  element  over  computations  that  employed  no 
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singular  element,  and  no  degradation  of  accuracy  was  observed  as  the  singular 
element  became  small.  This  study  leads  us  to  be  optimistic  that  the  small  sin¬ 
gular  element,  in  our  more  realistic  two-dimensional  cracks  between  different 
materials,  would  also  be  well  behaved.  To  confirm  this,  we  studied  a  model 
problem  on  two  different  meshes,  one  in  which  the  singular  element  is  the  same 
size  as  the  crack,  and  another  for  which  the  singular  element  is  20  times 
smaller.  In  addition,  we  have  calculated  the  stresses  using  an  entirely  dif¬ 
ferent  program  (FEAMOD)  that  employs  an  element  with  a  singular  transformation 
to  capture  the  square  root  behavior  (Stern,  1979)  (Reference  9).  The  differ¬ 
ence  between  this  element  and  the  singular  element  employed  in  FEAPICC  is  the 
form  of  the  basis  functions.  Recall  that  the  FEAPICC  singular  element  basis 
is  the  set  of  the  first  16  eigenfunctions  forming  the  exact  solution  to  the 
problem.  For  instance,  the  basis  includes  the  sinusoidal  log  r  behavior 
exhibited  by  cracks  between  different  elastic  materials. 

The  model  problem,  shown  schematically  in  Figure  5,  is  that  of  a  crack 
between  two  different  isotropic  materials  with  the  upper  material  100  times 
stiffer  than  the  lower  material.  Solutions  were  computed  with  FEAPICC  on  the 
coarse  grid  shown  in  Figure  6  and  on  the  fine  grid  in  Figure  7.  The  coarse 
grid  singular  element  is  shown  by  the  shaded  region  in  Figure  6.  The  fine 
grid  singular  elements  is  a  square  one-tenth  of  the  crack  length  on  edge.  The 
grid  for  the  FEAMOD  solution  is  shown  in  Figure  8  with  the  singular  element 
made  up  of  eight  triangular  elements  covering  the  area  inscribed  by  a  circle 
of  one  crack  length  in  diameter.  The  isoparametric  elements  in  the  FEAMOD 
solution  are  parabolic. 

Table  1  shows  the  stress  intensity  factors  for  each  of  the  three  numerical 
solutions  plus  the  exact  solution  of  the  problem  (Rice  and  Sih,  1965) 

(Reference  8).  The  two  FEAPICC  solutions  are  more  in  agreement  than  the  FEAMOD 
solution,  which  indicates  one  of  two  points.  This  difference  may  just  be 
indicative  that  the  stress 


Table  1.  Stress  Intensity  Model  Comparison 


I  Model  I 

i 

*»  I 

k2  1 

I  FEAPICC 

1 

(coarse)  | 

1.777  i 

-0.285  | 

j  FEAPICC 

(fine)  1 

1.809  1 

-0.291  | 

I  FEAMOD 

1 

1.5*2  | 

-0.300  j 

j  EXACT 

1 

1.793  | 

-0.262  | 
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Figure  5.  Crack  Between  Different  Isotropic  Materials 


Figure  6.  Coarse  FEAPICC  Grid 


Figure  7.  Fine  FEAPICC  Grid 


intensity  definitions  are  different  for  the  FEAMOD  and  FEAPICC  singular 
element,  or  it  may  show  a  real  difference  in  the  solutions.  Comparing  the  Y 
normal  stress  contours  in  the  vicinity  of  the  crack  on  the  exaggerated 
deformed  geometry  also  shows  some  differences  in  the  stress  field  suggesting 
that  the  FEAPICC  singular  element  gives  more  well-behaved  results.  Figures  9, 
10,  and  11  show  the  stress  contours  for  the  FEAPICC  coarse  and  fine  grid 
solutions  and  the  FEAMOD  solution,  respectively.  To  plot  the  results  for  the 
coarse  FEAPICC  and  FEAMOD  solutions,  the  stresses  and  displacements  on  the 
singular  element  have  been  interpolated  to  much  finer,  linear  elements.  The 
resulcs  of  the  fine  grid  singular  element  have  not  been  interpolated  to  yet  a 
finer  grid,  which  explains  the  weaker  contour  gradients  near  the  crack  tip. 

The  close  agreement  between  the  coarse  and  fine  grid  solutions  gives  us 
some  confidence  that  a  large  ratio  of  crack  length  to  element  length  will  not 
pose  any  problems.  This  may  be  especially  true  for  shock  wave  problems  in 
which  the  shock  width  may  be  much  smaller  than  the  delamination  length. 

Hence,  the  stress  concentration  effects  are  localized  near  the  crack  tip, 
since,  to  the  shock  wave,  most  of  the  delamination  looks  like  one  uninterrupted 
free  surface.  Indeed,  numerical  experiments,  to  be  described  later,  indicate 
such  a  localization. 

3 . 2  Thermal  Stress  Validation 

To  test  the  implementation  of  the  thermal  stresses  in  the  FEAPICC  code,  the 
singular  finite  element  solution  is  compared  to  a  simple  analytical  solution 
that  has  no  singularity.  The  singular  solution  should  be  in  close  agreement 
except  near  the  crack.  In  the  problem  shown  in  Figure  12,  the  temperature 
varies  linearly  across  the  plate  in  plane  stress.  With  the  loading  as  shown, 
the  closed  form  solution  for  an  isotropic,  homogeneous  plate  is 

u  =  -VO  x/E  +  a  T  xy/L 

° 2  2  (33) 

v  =  a  y/E  +  a  TQ(y  -  x  )/2L 

The  finite  element  solution  has  a  small  crack  of  length  L/20  starting  at 
y  =  L/2,  x  =  0.  Except  in  the  vicinity  of  this  crack,  the  solution  compares 
closely  to  the  analytical  solution.  Table  2  compares  the  two  solutions  for 
various  points  on  the  plate.  The  displacement  at  the  crack  is  an  average  of 
the  top  and  bottom  surfaces  of  the  crack. 


Figure  9.  Normal  Stress  FEAPICC  Coerse  Grid 
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Figure  10.-  Normal  Stress  FEAPICC  Fine  Grid 
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Table  2.  Static  Thermoelaotic  Teat  Case  Comparison 

(L  =  20,  a  =  1.0,  E  =  1.0,  \>  =  0.25,  AT  *  1.0) 
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1  u 
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1 

1 

10 

1 

0 

1  -2.50  | 

-2.41 

1 
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1 
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1 

1 

0 

1 

5 

1  o  I 

0 

1 
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5.57 
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5 
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1 

1 
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5 
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1 

3.12 

1 
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1 

1 

0 

1 

10 
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0 

1 
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1 

12.89 

1 

1 

5 

1 

10 
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0.94 

1 
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1 
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1 

1 

10 

1 

10 
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1 
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1 
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1 

1 

0 

1 

15 
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1 
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1 
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2.42 
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1 

20.26 

1 

1 

10 

1 
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1 
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1 
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0 

1 

20 
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1 

5 

1 
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1 
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1 

20 
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1 

27.5 

1 
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3.3  Static  Composite  Plate  Validation  Test 

In  this  section,  we  investigate  the  static  response  of  a  multi-ply  plate 
with  a  centered  delamination  that  has  dimensions  much  larger  than  the  ply 
thickness.  The  objective  is  to  test  the  proper  implementation  of  the  planar 
stress/strain  relations  with  a  more  realistic  material  model.  The  grid  and 
loading  for  these  tests  is  shown  in  Figure  13.  The  plate  consists  of  eight 
plies  with  two  elements  within  each  ply.  The  material  properties,  typical  of 
graphite-epoxy,  are  given  in  Table  3.  Two  composite  plates  are  considered:  a 
center  symmetric  plate  (s)  with  fibers  lying  in  the  successive  plies  starting 
from  the  bottom  at  0  =  0°,  45° ,  90°,  -45°,  -45°,  90°,  45°,  and  0°  and  an 
unsymmetric  plate  (u)  with  successive  plies  at  0=0°,  45°,  90°,  -45°,  0°, 

45°,  90°,  -45°.  The  two  different  composites  have  been  subjected  to  three 
different  uniform  temperature  distributions: 

(1)  plate  at  reference  temperature,  AT  =  0  (off); 

(2)  plate  below  reference  temperature,  AT  =  -0.1°F  (cold);  and 

(3)  plate  above  reference  temperature,  AX  =  0.1°F  (hot), 
and  two  pressure  loads: 

(1)  no  mechanical  load  (off)  and 

(2)  a  load  of  1  psi  (on). 
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L  -  0.25" 
a  -  0.05" 
b  -  0.04" 
t  -  0.005" 
p  ■  1  psi 
d  -  0.025 


Figure  13.  Composite  Grid 


Table  3.  Material  Properties  Along  the  Principal  Coordinates 
of  the  Fiber  Direction 


Ei  -  1.8xl07  psi 
\>12  =  Vi3  =  0.34 
G12  =  G13  =  0.95x10^  psi 
ai  =  2.0xl0“7  (“F)”1 


e2  “  E3  *  1.4  xlO^  psi 
\>23  ”  0.4  p  *  0.055  lb/in^ 
G23  *  0.5x10^  psi 
02  =  a3  =  1.6x10“5  (°F)-1 


The  magnitudes  of  the  thermal  loads  are  chosen  so  that  they  are  approxi¬ 
mately  the  same  magnitude  as  the  pressure  loads.  Table  4  shows  the  stress 
intensity  factors  for  each  ply  stacking,  pressure,  and  thermal  loading. 

It  is  easy  to  show  that  the  stress  intensity  factors  behave  as  a  linear 
function  of  the  mechanical  load  and  thermal  load.  Also,  the  strain  energy 
release  rate  is  a  quadratic  function  of  the  mechanical  load  and  the  thermal 
load.  The  numerical  results  in  Table  4  indicate  that  the  strain  energy 

I  release  rate  shows  a  strong  coupling  between  the  mechanical  load  and  the 

thermal  load.  This  implies  that  the  strain  energy  resulting  from  a  combined 
mechanical  and  thermal  load  is  not  simply  the  sum  of  the  strain  energies  of 
each  load  independently. 
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Table  4.  Stacie  Composite  Ply  Stress  Intensity  Pactors 
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1 
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1 
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3.4  Dynamic  Composite  Plate  Validation  Test 

Before  moving  to  the  actual  simulation  of  the  stress  pulses  generated  by 
laser  heating,  we  examine  in  this  section  the  stress  wave  propagation  generated 
by  a  pulsed  tensile  loading  on  the  surface  of  the  plate.  We  chose  a  pulse 
loading  that  is  characteristic  of  that  produced  by  pulsed  laser  after  the  back 
free  surface  reflection.  The  grid  shown  in  Figure  13  and  the  symmetric  stack¬ 
ing  sequence  and  material  properties  described  in  Section  3.2  are  used  in  this 
simulation.  With  these  material  properties  we  find  a  maximum  longitudinal  wave 
speed  of  3.6x10^  in/s  along  the  fiber  direction.  The  transverse  longitudinal 
wave  speed  is  about  1.0x10^  in/s,  while  the  shear  wave  speed  ranges  from 

0.6x10^  ir./s  in  the  0°  ply  to  0.8x10^  in/ s  in  the  90°  ply.  Since  the 

-3  .  .  -9 

smallest  grid  dimension  is  2.5x10  in,  a  step  size  of  5.0x10  s  allows 
the  fastest  wave  to  traverse  less  than  one  element  in  one  step. 

In  our  simulations  we  noted  severe  oscillation  problems  with  a  square 
tensile  loading  pulse  wave.  These  oscillations  are  characteristic  of  the 
element's  inability  to  resolve  the  sharp  front  of  the  wave.  To  remove  this 
effect,  we  selected  a  hat-shaped  pulse  loading  in  which  the  loading  increases 
linearly  with  time  until  some  maximum  is  reached  and  then  falls  linearly  back 
to  zero.  The  time  to  reach  maximum  loading  is  set  to  2.0x10  ^  s.  This  produces 


a  stress  wave  within  the  body  that  rises  and  falls  within  eight  elements  or 
four  ply  widths,  and  the  oscillations  are  only  modest.  Ratios  of  half¬ 
wavelength  to  element  size  smaller  than  four  show  marked  oscillations 
regardless  of  the  step  size  using  the  Newmark  method.  The  Wilson  method  with 
its  inherent  numerical  damping  proved  to  be  only  a  slightly  better  option, 
since  little  damping  was  apparent  fcr  small  time  step  sizes  and  too  much 
damping  or  loss  in  accuracy  was  observed  for  the  larger  step  sizes.  These 
wave  resolving  restrictions  have  some  important  implications  especially  for 
sharp  laser  pulses  in  three-dimensional  simulations,  which  we  discuss  further 
in  Section  6 . 

In  Figures  14  through  23,  we  show  the  stress  contours  at  different  times 
for  the  simulation  of  the  symmetric  plate.  These  contours  are  plotted  on  the 
exaggerated  deformed  geometry.  Since  the  tensil?  load  at  a  maximum  of  1.0x10^ 
psi  is  applied  uniformly  across  the  top  surface,  the  stress  wave  is  very 
uniform  before  it  interacts  with  the  delamination.  In  Figure  14,  we  see  the 
normal  transverse  wave  at  a  time  just  before  its  head  reaches  the  delamination 
and  before  the  end  of  the  boundary  loading  pulse.  Note  that  the  delamination 
(marked  by  the  line  ending  with  the  asterisk)  is  assumed  closed  initially. 
Although  the  delamination  is  initially  in  the  center,  the  tensile  deformation 
in  the  wave  has  stretched  the  top  plies  in  this  exaggerated  plot.  In 
Figure  15,  the  transverse  normal  stress  is  shown  as  the  maximum  stress  passes 
in  front  of  the  crack.  We  see  that  the  delamination  has  opened  almost 
uniformly  across  its  length.  Also,  we  see  the  first  signs  of  a  compression 
wave  reflecting  off  the  free  surface  of  the  delamination.  Five  hundredths  of  a 
microsecond  later  (Figure  15)  the  compression  wave  is  more  fully  developed,  and 
still  later  (Figure  17)  we  see  the  interaction  of  the  compression  and  tensile 
waves  with  the  upper  and  lower  surfaces.  The  stress  intensity  factors  labeled 
on  each  plot  show  that  these  peak  wher  both  large  compression  and  tensile 
stresses  are  near  the  crack  tip.  The  shear  wave  propagates  in  a  characteristic 
antisymmetric  profile  as  shown  in  Figures  18  through  20.  The  time  sequence  of 
the  longitudinal  normal  stress  shown  in  Figures  21  through  23  indicates  the 
dramatic  effects  of  fiber  stiffening.  In  the  90°  plies  the  fibers  are  running 
parallel  to  the  long  axis.  Since  the  waves  travel  over  3  times  faster  in  90° 
plies,  we  see  that  the  contours  seem  to  spurt  out  in  these  plies. 
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Figure  14.  Normal  Stress  (t  -  0.15  ys) 


LOADCASE:  SO  TIMESTEP:  SO  TINE:  0.0000003 
FRAME  OF  REF:  6L0SAL 

STRESS  -  Y  NIM:  -2.26E+04  MAX:  1.00E409 


Kx  -  8. A  x  10J 


Kj  •  -3.1  x  10* 


Figure  15.  Normal  Stress  (t  ■  0.3  ys) 
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Figure  16.  Normal  Stress  (t  “  0.35  ye) 
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Figure  17.  Normal  Stress  (t  *  0.4  ys) 
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Figure  20.  Shear  Stress  (t  ■  0.4  us) 


Figure  21.  Longitudinal  Normal  Stress  (t  «  0.3  ys) 
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A.  SIMULATION  OF  THE  STRESS  WAVE  INITIATION  AND  PROPAGATION  IN  ONE  DIMENSION 

A  number  of  finite  difference  solutions  to  the  complete  problem  of  energy 
deposition  and  subsequent  motions  were  obtained.  The  code  WONDY  (Lawrence  and 
Masur,  1971)  (Reference  3)  was  utilized  to  obtain  these  results.  A  20-layer 
composite  laminate  was  divided  into  a  total  of  400  zones,  with  20  zones  in 
each  layer.  The  total  thickness  was  0.1  inch,  or  0.254  cm;  each  layer  was 
then  one-tenth  of  that  thickness.  The  mass  density  of  the  material  was  1.568 

3 

gm/cm  ,  which  is  typical  of  an  epoxy-graphite  material  (Froula  et  al.,  1980) 
(Reference  2).  The  sound  speed  was  3.02x10"*  cm/s  (Lee,  1979)  (Reference  4). 

A  vapor  transition  was  included;  the  vaporization  energy  was  set  at  10^ 
ergs/gm  (Lee,  1979)  (Reference  4). 

For  the  solid  region,  a  Mie-Gruneisen  equation  of  state  was  used.  The 

reference  Gruneisen  parameter  was  equal  to  0.3  (Froula  et  al.,  1980) 

(Reference  2).  A  linear  wave  speed  versus  particle  velocity  form  was  adopted, 

with  a  slope  of  0.65  (Lee,  1979)  (Reference  4).  Thus,  the  material  did  not 

have  a  constant  wave  speed,  but  the  model  had  nonlinearities. 

For  the  stress  deviator  response,  a  linear  elastic  model  with  a  von  Mises 

plastic  yield  was  used.  The  Poisson's  ratio  was  0.44  (Lee,  1979) 

(Reference  4),  which  gives  a  unaxial-strain  wave  speed  of  3.26x10"*  cm/s. 

Since  the  behavior  for  the  initial  energy  deposition  and  subsequent  wave 

motions  is  dominated  by  the  hydrodynamic  model,  it  was  not  necessary  to  use  an 

anisotropic  model  for  this  part  cf  the  problem. 

9  2 

The  yield  stress  was  0.5x10  dynes/cm  .  Failure  in  tensile  spall  was 
9  2 

assumed  to  occur  at  10  dynes/cm  (14,500  psi). 

With  a  longitudinal  speed  of  3, 26x10"*  cm/ s ,  and  a  total  thickness  of 
0.254  cm,  the  total  plate  transit  time  is  0.78x10  ^  s ,  or  just  under  a  micro¬ 
second.  The  energy  deposition  time  was  set  at  one-fortieth  of  a  microsecond, 
or  0.25x10  ®  s.  The  energy  was  deposited  at  a  constant  rate  over  this  time. 
For  the  spatial  distribution,  the  linear  shape  discussed  in  the  previous 
section  was  used  as  an  approximation  to  the  exponential  shape  characteristic 
of  actual  depositions.  Thus,  the  maximum  energy  density  occurred  at  the  front 
surface,  and  that  energy  density  decreased  linearly  to  zero  at  a  total  depth 

of  one-twentieth  of  the  total  thickness  (0.0127  cm  or  0.005  inch).  The  total 
8  2 

flux  was  2x10  ergs/cm  .  With  the  deposition  thickness  of  0.0127  cm  and  the 
mass  density  of  1.568  gm,  this  gives  an  average  initial  energy  density  of 
10^  ergs/gm.  At  the  surface,  the  energy  density  is  therefore  2x10^  ergs/gm. 
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Since  the  Griineisen  is  0.3,  an  instantaneous  deposition  of  this  energy  would 

10  9  2 

create  an  initial  stress  of  P  =  (1 . 568) ( 0 . 3)(2xl0  )  =  9.5x10  dynes/cm  ,  or 

just  under  10  kilobars,  with  a  pulse  width  of  0.0127  cm.  In  fact,  because  of 

the  finite  deposition  time,  the  actual  pulse  width  is  slightly  higher,  and  the 

initial  stress  about  half  that  value. 

As  discussed  in  Section  1.3,  an  initial  triangular  compressive  pulse  will 

generate  a  tensile  following  "tail"  as  it  moves  away  from  the  free  surface. 

As  that  tail  builds  in  magnitude,  tensile  spall  of  the  front  surface  will 

occur.  Additional  spall  will  occur  until  the  final  tensile  magnitude  decreased 

to  under  the  spall  strength  of  1  kilobar. 

Figure  24  shows  the  actual  resultant  wave  at  the  time  of  10  ^  s,  well 

after  the  final  wave  has  moved  away  from  the  free  front  surface.  The  main 

3 

compressive  pulse  is  about  3  kilobars  (40x10  psi)  at  this  time  and  is 
indeed  almost  triangular  in  shape.  At  times  immediately  after  the  energy 
deposition,  parts  of  the  front  surface  did  spall  as  the  tensile  wave  began  to 
grow.  Ultimately,  a  total  thickness  one-half  the  deposition  depth  was  spalled. 
The  tensile  following  wave,  which  was  almost  a  kilobar  at  the  time  of  final 
spall,  has  at  this  time  decreased  to  about  two-thirds  of  a  kilobar. 

This  result  illustrates,  the  important  effects  of  the  material 
nonlinearities  and  especially  the  effects  of  the  spall  response.  It  is  u.  j 
worthwhile  noting  that,  in  this  particular  problem,  the  initial  energy 
densities  were  below  the  vaporization  energy  (10^  ergs/gm)  by  a  factor  of  5. 

If  the  energy  flux  had  been  greater,  the  material  near  the  free  surface  would 
have  been  vaporized,  and  the  tensile  strength  would  have  decreased  to  zero. 

The  wave  shown  in  Figure  24  propagates  into  the  interior  of  the  plate. 

Since  the  material  is  nonlinear,  the  stress  pulse  will  decay  in  magnitude  and 
broaden  in  width.  For  example,  Figure  25  shows  the  resultant  wave  at  the  time 
of  0.7x10  ^  s,  just  before  it  reaches  the  back  free  surface  of  the  plate. 

At  this  time,  the  peak  compressive  wave  has  decayed  to  30,000  psi,  and  the 
tensile  portion  has  decayed  to  less  than  6000  psi.  The  pulse  is  almost  twice 
as  wide  as  it  was  initially. 
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Figure  24.  Initial  Fully  Developed  Stress  Wave  (t  ■  0.1  \m) 


Figure  25.  Stress  Have  Just  Before  Free  Surface  Reflection  (t  ■  0.7  ye) 
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5 .  SIMULATION  OF  THE  STRESS  WAVE  INTERACTION  WITH  A  DELAMINATION 

As  the  stress  wave  reflects  from  the  wall  one  would  expect  the  surface 
plies  to  spall.  To  simulate  such  an  event,  we  model  the  four  plies  adjacent  to 
the  surface.  In  the  middle  of  these  four  plies  is  a  small  delamination  with  a 
half-length  of  0.025  inch.  The  grid  is  the  same  as  that  used  in  the  previous 
composite  simulations  (Figure  13),  but  the  dimensions  are  reduced  by  half  to 
give  four  elements  per  ply.  The  initial  conditions  are  taken  from  the  one¬ 
dimensional  simulation  at  a  time  just  before  the  reflection  of  the  wave  but 
just  after  the  compressive  part  of  the  wave  passes  the  delamination.  The  top 
and  bottom  surfaces  are  both  considered  stress-free.  Since  the  stress  wave 
lies  ithin  about  two  plies,  the  bottom  surface  observes  essentially  no  stress 
from  the  time  the  compressive  wave  passes  over  until  the  time  the  reflected 
tensile  wave  returns.  It  is  during  this  time  frame  that  we  have  conducted  a 
series  of  simulations  varying  the  ply  stacking  sequence,  plate  temperature, 
crack  geometry,  and  boundary  restraints.  Table  5  gives  a  summary  of  the  test 
cases  and  the  maximum  fracture  intensity  values  observed  over  each  simulation. 
In  the  first  six  cases  we  vary  the  stacking  sequence  and  observe  that  the 
maximum  fracture  parameters  vary  little  from  case  to  case.  Even  though  the 
loading  is  not  symmetric,  the  plates  with  delaminations  between  similiar  plies 
show  no  mode  II  effects.  Pistes  with  delaminations  between  dissimilar  plies 
show  a  small  mode  II  effect,  but  still  it  is  minor  compared  to  a  mode  I 
fracture.  The  mode  II  strain  energy  release  rate  is  especially  small  because 
the  differential  in  the  deformation  tangent  to  the  crack  is  small.  Since  the 
mode  I  strain  energy  release  rate  is  defined  primarily  by  the  Young's  modulus 
of  the  matrix,  this  parameter  varies  little  with  the  ply  sequence. 

Wang  et  al.  (1980)  (Reference  12)  found  that,  since  the  curing  temperatures 
are  so  high  (300°F),  the  thermal  stress  at  room  temperature  is  an  important 
consideration  when  evaluating  the  strain  energy  release  rate.  To  test  this 
effect,  in  case  7,  the  plate  is  assumed  to  be  200°F  colder  than  the  reference 
temperature.  In  comparison  to  case  1,  a  significant  change  in  the  strain 
energy  release  rate  is  seen,  indicating  that  the  ambient  temperature  of  the 
plate  plays  an  important  role  in  defining  its  fracture  resistance.  In  this 
simulation,  we  assumed  that  the  temperature  of  the  plate  suddenly  dropped 
200#F.  This  may  give  rise  to  some  unwanted  dynamical  effects  near  the 
constrained  regions  of  the  plate.  Since  the  delamination  is  removed  from 
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Table  5.  Dynamic  Parametric  Tests 
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these  constraints,  we  expect  that  these  effects  are  minor.  In  case  8,  the 
delamination  is  assumed  to  reside  on  a  free  edge  of  the  plate.  Note  that  the 
fracture  parameters  changed  very  little.  This  is  somewhat  surprising  since, 
in  static  tests,  the  edge  delamination  is  much  more  apt  to  fracture. 

The  reason  these  dynamic  tests  do  not  show  this  behavior  can  again  be 
attributed  to  the  pulse  wavelength  to  delamination  length  ratio.  The  tip  of 
the  delamination  is  initially  unaware  that  the  free  surface  is  nearby.  The 
free  edge  effects  would  be  apparent  only  for  delaminations  on  the  order  of  the 
wavelength  of  the  pulse. 

In  cases  9  through  11,  the  boundary  is  constrained  so  that  only  that 
portion  of  the  boundary  above  the  delamination  is  free.  As  the  wave  reflects 
off  this  back  surface,  a  shear  wave  will  ensue  from  the  region  near  the 
constrained/unconstrained  boundary.  This  shear  wave  propagates  into  the 
interior  and  interacts  with  the  delamination  to  produce  excessive  stresses. 
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Since  the  shear  wave  is  partially  transmitted  via  the  fibers,  one  would  expect 
to  see  variations  in  the  fracture  parameters  depending  on  the  stacking 
sequence.  The  tests  conducted  verify  this  sensitivity.  Since  the  shear  waves 
travel  alowcr  than  the  normal  stress  waves,  the  maximum  value  of  the  mode  II 
fracture  parameters  peak  at  a  later  time  than  the  mode  I  parameters. 

In  all  the  simulations,  the  maximum  fracture  intensity  conditions  occurred 
after  the  maximum  of  the  tensile  stress  wave  had  passed  slightly  below  the 
crack  and  the  compressive  wave  had  partially  reflected  from  the  crack  f*ee 
surface.  The  critical  strain  energy  release  rate  for  unstable  growth  is  still 
in  question,  but  under  static  loading  tests  Wang  et  al.  (1980)  (Reference  12) 
quoted  a  value  of  0.8  lb/in.  This  would  imply  that  the  spalling  occurs  well 
before  the  maximum  intensity  conditions  at  a  time  just  before  the  maximum 
stress  reaches  the  crack.  In  such  a  scenario,  most  of  the  energy  would  remain 
as  kinetic  energy  in  a  spall  the  size  of  which  would  be  related  to  the  stress 
wavelength.  This  would  effectively  damp  the  wave  in  the  main  body  by 
preventing  this  energy  from  being  converted  to  strain  energy.  The  spall  would 
travel  at  a  high  rate  of  speed  away  from  the  main  body.  The  velocity  vector 
plot  in  Figure  26  shows  speeds  on  the  order  of  150  mph  at  the  time  of  maximum 
fracture  intensity. 

None  of  the  simulations  indicate  that  the  initial  tensile  tail  plays  much 
of  a  role  in  the  fracture  of  the  composite.  This  may  not  be  true  for  loading 
conditions  other  than  the  magnitude  and  profile  assumed  in  the  one-dimensional 
simulations  in  Section  4. 

Plots  of  the  transverse  normal  and  shear  stress  contours  within  the 
singular  element  (0.0025  inch  on  edge)  given  in  Figures  27  and  28  show  a  very 
typical  pattern.  These  plots  are  at  a  time  of  maximum  fracture  intensity.  As 
the  wave  interacts  with  the  delamination,  the  region  near  the  crack  tip  could 
behave  plastically.  With  this  in  mind,  the  contours  of  von  Hises  stress 
plotted  in  Figure  29  give  an  indication  of  the  size  of  the  plastic  region. 
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Figure  26.  Velocity  Vector  at  Maximum  Stress  Intensity  (t  -  0.16  ys) 


Figure  27.  Normal  Stress  in  the  Singular  Element  (t  ■  0.16  ye) 
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6.  UNIFIED  MODEL  FEASIBILITY  STUDY 

The  long-range  goal  of  this  study  is  the  development  of  a  complete  finite 
element  program  for  rapid  thermal  loading  that  is  capable  of  an  accurate  and 
complete  thermodynamic  description  similar  to  the  finite  difference  code  used 
in  this  study.  As  an  interim  measure  we  have  used  the  finite  difference  model 
to  describe  the  energy  deposition  aspects  of  the  problem  and  to  follow  the 
resulting  waves  until  they  decay  into  the  lower  temperature  solid  behavior. 

At  this  level,  a  transfer  to  the  finite  element  code  is  made  to  study  the  sub¬ 
sequent  interaction  with  a  pre-existing  delamination.  Such  a  treatment  cannot, 
however,  properly  treat  the  coupling  of  the  effects  of  the  crack  on  the  stress 
wave  generation,  which  could  be  important  for  delaminations  near  the  irradiated 
surface.  In  addition  to  these  coupling  effects,  each  of  the  individual  models 
used  in  this  Phase  I  study  are  inadequate,  in  certain  respects,  for  describing 
the  high-energy,  high-stress  behavior  of  the  composite  material.  This  is  in 
part  due  to  the  lack  of  experimental  data  related  to  composites  in  this  regime. 
Nevertheless,  certain  aspects  could  be  improved  with  current  knowledge.  For 
example,  in  the  current  finite  element  model,  the  composite  is  treated  as  a 
linearly  elastic  anisotropic  material.  Certainly,  as  the  plate  is  heated,  an 
epoxy  matrix  would  undergo  some  phase  transition  that  alters  its  elastic 
properties,  while  the  graphite  fibers  might  retain  their  elastic  properties  at 
the  same  temperature.  The  result  would  be  some  anisotropic  medium  that  behaves 
like  a  fluid  in  certain  directions.  Similarly,  after  the  wave  has  moved  out 
of  the  high-temperature  region,  the  high  stresses  might  cause  the  matrix  to 
deform  plastically. 

Certain  algorithmic  features  need  to  be  employed  to  allow  for  more  than 
one  delamination,  create  new  delaminations,  and  allow  for  partial  or  full 
closing  of  delaminations  depending  on  the  nature  of  the  stress  wave.  The  logic 
for  closing  delaminations  may  be  conveniently  expressed  in  the  finite  element 
context  via  Lagrangian  multipliers.  These  multipliers  represent  the  contact 
force  between  the  sides  of  the  delamination.  In  addition  to  interlaminar 
delaminations,  the  model  should  consider  transverse  cracks,  which  may  be 
important  when  normal  operating  loads  are  also  considered  in  the  analysis. 

In  the  current  work,  the  stress  wave  generation  was  carried  out  in  one 
dimension.  This  limits  the  analysis  to  the  generation  of  longitudinal  waves. 
The  next  step  would  be  to  carry  out  the  analysis  in  two  dimensions  under  plane 
stress/strain  or  axisymmetric  assumptions.  Strictly  speaking,  the  anisotropic 


and  inhomogeneous  properties  of  the  composite  plate  rule  out  any  assumptions 
of  axisymmetry  except  in  idealized  cases.  Similarly,  plane  stress/strain 
assumptions  are  applicable  only  by  restricting  the  nature  of  the  laser  area 
deposition.  The  problem  is  truly  three-dimensional,  but  the  merit  of  three- 
dimensional  calculations  must  be  weighed  with  their  complexity  and  the  state 
of  the  art  of  the  physical  model.  Three-dimensional  simulations  may  be 
feasible  under  certain  settings.  Since  the  grid  size  must  be  related  to  the 
wavelength  of  the  stress  pulse,  which  is  given  by  the  laser  pulse  duration, 
short  pulses  imply  fine  grids.  For  the  same  amount  of  work,  one  might  model 
long  pulses  over  a  large  area,  or  short  pulses  over  a  small  area. 

Given  that  the  physical  models  describing  composites  under  extreme  loading 
conditions  are  still  in  their  infancy,  it  is  absolutely  imperative  that  the 
computer  model  be  modular  and  structured.  This  will  enable  the  model  to 
evolve  as  experimental  evidence  and  physical  models  become  available.  All  the 
essential  processes  need  to  be  at  least  identified  early  in  the  model 
development  to  form  a  framework  to  build  on. 

To  this  end,  the  finite  element  method  offers  just  such  a  modular  struc¬ 
ture.  Codes  such  as  DYNA2D  and  DYNA3D,  available  in  the  public  domain,  are 
based  on  a  finite  element  formulation,  which  includes  finite  deformation, 
large  strain,  and  critical  bulk  viscosity  to  control  shock  waves.  The  codes 
include  a  mesh  rezoner  to  maintain  a  proper  mapping  in  these  Langrangian 
calculations . 

Unlike  the  FEAPICC  codes,  the  formulation  is  explicit.  This  implies  that 
no  matrices  need  to  be  assembled  or  eliminated,  reducing  the  computer  time  and 
memory.  Since,  for  time  accuracy,  the  step  sizes  must  be  set  below  the 
stability  limit,  the  explicit  formulation  offers  no  drawbacks  in  this  regard. 
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7 .  CONCLUSIONS 

To  summarize,  we  have  addressed  the  problem  of  rapid  deposition  of  thermal 
energy  on  the  surface  of  a  composite  plate.  Several  general  aspects  of  the 
problem  have  been  described,  and  the  assumptions  regarding  the  analysis  have 
been  enumerated.  Existing  codes  have  been  modified  to  solve  the  stress  pulse 
generation  and  its  subsequent  interaction  with  a  delamination  within  the 
composite  plate.  The  finite  difference  model  for  simulation  of  the  stress 
wave  generation  includes  the  thermomechanical  coupling  and  a  general  equation 
of  state  model  to  account  for  phase  changes.  The  finite  element  model  accounts 
for  the  stress  singularity  of  a  delamination  within  the  anisotropic  composite 
plies.  This  model  has  been  modified  to  include  thermal  stresses,  to  calculate 
the  strain  energy  release  rates,  and  to  account  properly  for  the  plane  strain 
assumptions.  Tests  have  been  conducted  to  ensure  that  the  modifications  have 
been  implemented  correctly  and  that  the  model  is  applicable  to  the  problem  we 
address . 

Our  analysis  shows  that  the  strain  energy  is  a  quadratic  function  of  the 
mechanical  and  thermal  loading  with  a  strong  coupling  between  the  mechanical 
and  thermal  parts.  While  the  study  of  Wang  et  al.  (1980)  (Reference  12) 
implies  the  contrary,  the  strain  energy  resulting  from  a  combined  mechanical 
and  thermal  load  is  not  generally  the  sum  of  the  strain  energies  of  each  load 
independently . 

The  results  of  the  analysis  of  the  stress  wave  generation  indicate  that 
the  ensuing  wave  will  have  both  compressive  and  tensile  parts.  The  magnitude 
of  the  tensile  part  depends  on  the  tensile  strength  or  fracture  toughness  of 
the  surface  plies.  As  the  wave  propagates  into  the  interior  cf  the  plate,  the 
nonlinearities  damp  and  spread  this  wave. 

We  have  modeled,  using  the  singular  finite  element  method,  the  spalling 
event  as  the  stress  wave  reflects  from  the  back  surface  and  interacts  with  an 
existing  delamination.  Our  analysis  indicates  that  the  fracture  parameters  are 
insensitive  to  the  stacking  sequence  when  the  delaraination  interacts  with  a 
pure  tension  wave.  This  is  because  the  stress  is  transmitted  primarily  through 
the  matrix.  If  a  shear  wave  results  from  the  loading  or  boundary  constraints, 
then  the  stacking  sequence  has  an  appreciable  effect  on  the  fracture  param¬ 
eters,  since  the  shear  wave  is  transmitted  via  the  fiber/matrix  combination. 

The  ambient  temperature  of  the  plate  also  proved  to  be  an  important  considera¬ 
tion  in  the  fracture,  since  the  thermal  deformation  resulting  from  the  curing 
of  the  composite  can  be  significant. 
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There  are  many  aspects  of  this  model  that  need  more  attention.  More 
experimental  studies  and  models  are  needed  to  characterize  the  behavior  of 
composites  under  extreme  loads.  We  recommend  that  the  models  for  stress 
generation  and  for  its  subsequent  interaction  with  a  delamination  be  unified 
under  one  finite  element  framework.  The  inherent  modular  structure  of  the  FEM 
will  lend  itself  to  modifications  as  experimental  data  become  available.  Such 
a  code  could  also  help  guide  the  experimental  program  and  aid  the  experimen¬ 
talists  in  interpreting  their  results. 
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